#setwd("/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/1000samples")
#setwd("/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22")
setwd("/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22")
map<-read.table("map.txt",sep="\t",head=T)
segdat<-read.table("segdat.txt",sep="\t",head=T)
outdir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/result_new_paper_graph_11_01_2011"
setwd(outdir)
#png(file="chr22_median_dosage_diff_1000samples.png",res=100, width=1000, height=800)
#png(file="chr22_median_dosage_diff_allsamples.png",res=100, width=1000, height=800)
png(file="chr22_median_dosage_diff_allsamples.png",res=100, width=1600, height=1600)
yrange<-c(0,0.20)
plot(lowess(map$bp, map$meddiff.seed,f=0.00001), type='l', ylim=yrange, lty=1,main="Median of Absolute Value of Dosage Difference", xlab="BP", ylab="Median Absolute Difference")
lines(lowess(map$bp, map$meddiff.4mb,f=0.00001), lty=2, col="brown")
lines(lowess(map$bp, map$meddiff.2mb,f=0.00001), lty=3, col="blue")
lines(lowess(map$bp, map$meddiff.1mb,f=0.00001),lty=4, col="red")
lines(lowess(map$bp, map$meddiff.halfmb,f=0.00001),lty=5, col="green")
lines(lowess(map$bp, map$meddiff.quatermb,f=0.00001), lty=6, col="cyan")
lines(lowess(map$bp, map$meddiff.halfquatermb,f=0.00001), lty=7, col="purple")
lines(lowess(map$bp, map$meddiff.quaterquatermb,f=0.00001), lty=8, col="magenta")
#names(segdat)<-c("break.","seg.quaterquatermb","seg.halfquatermb","seg.quatermb","seg.halfmb","seg.1mb","seg.2mb","seg.4mb")
lines(lowess(map$bp[!is.na(map$r2)], yrange[2]-yrange[2]*map$r2[!is.na(map$r2)],f=0.00001), lty=9, col="orange")
abline(v=segdat$seg.4mb, col="brown", lty=2)
abline(v=segdat$seg.2mb, col="blue", lty=3)
abline(v=segdat$seg.1mb, col="red", lty=4)
abline(v=segdat$seg.halfmb, col="green", lty=5)
abline(v=segdat$seg.quatermb, col="cyan", lty=6)
abline(v=segdat$seg.halfquatermb, col="purple", lty=7)
abline(v=segdat$seg.quaterquatermb, col="magenta", lty=8)
legend("top", col=c("orange","black", "brown","blue","red", "green","cyan","purple","magenta"), lty=c(9,1,2,3,4,5,6,7,8,9), c("Scaled 1-r2,\nPeaks are\nareas of low quality","Different\nSeed", "4MB", "2MB", "1MB","0.5MB","0.25MB","0.125MB","62.5KB"), horiz=T,cex=0.75)
dev.off()
